Introducción al Análisis de la Microbiota Vaginal con Tidytacos
Monica R. Ticlla
2023-12-11
Esta guía explica cómo realizar algunas tareas comunes de análisis de datos de recuento de secuenciación microbiana utilizando tidytacos. Ilustraremos utilizando un conjunto de datos con muestras de microbioma humano del tracto vaginal, tomadas de este artículo por Lebeer et al.
Tidytacos: la filosofía del paquete
Un objeto tidytacos es simplemente una lista de tres tablas:
counts: estos son los recuentos de lecturas para cada taxón (OTU/ASV/filotipo) en cada muestra. Cada fila representa dicho recuento de lecturas.
samples: esta tabla contiene los metadatos de muestra. Cada fila representa una muestra.
taxa: esta tabla contiene la taxonomía y otros metadatos de los taxones. Cada fila representa un taxón.
El paquete se llama tidytacos porque cada una de las tablas está ordenada (en inglés tidy): cada fila representa una observación y cada columna una variable (puedes encontrar más información sobre la ordenación de datos en este enlace).
Configuración
En caso de que aún no hayas instalado tidytacos, puedes instalarlo usando devtools:
install.packages("devtools")
devtools::install_github("LebeerLab/tidytacos")Para esta guía, sólo necesitamos cargar tres paquetes: tidytacos (por supuesto), el conjunto de paquetes tidyverse, y VTutorials.
library(tidyverse)
library(tidytacos)
library(VTutorials)Exploremos el dataset vdata
Nuestro conjunto de datos de ejemplo está disponible en el paquete VTutorials y no es necesario importarlo ni convertirlo. Nuestro set de datos de ejemplo Se llama “vdata”.
“vdata” es un objeto de tipo tidytacos y contiene 3 tablas: “counts”, “samples”, y “taxa”.
summary(vdata)## Length Class Mode
## counts 3 tbl_df list
## samples 28 tbl_df list
## taxa 12 tbl_df listComenzamos inspeccionando la tabla de muestras:
glimpse(vdata$samples)## Rows: 200
## Columns: 28
## $ sample <chr> "SAMEA13411869", "SAMEA13410951",…
## $ sample_id <chr> "s16", "s100", "s114", "s170", "s…
## $ Technical.run_20201009_isala_cross_1 <lgl> FALSE, TRUE, FALSE, FALSE, TRUE, …
## $ Technical.run_20201016_isala_cross_2 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201023_isala_cross_3 <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, …
## $ Technical.run_20201030_isala_cross_4 <lgl> FALSE, FALSE, FALSE, TRUE, FALSE,…
## $ Technical.run_20201106_isala_cross_5 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201120_isala_cross_6 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201204_isala_cross_8 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201211_isala_cross_9 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20210119_isala_cross_7 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.libsize <dbl> 37957, 22782, 32150, 7707, 24123,…
## $ General.Age <dbl> 24, 60, 24, 33, 29, 24, 23, 42, 3…
## $ Health.BMI <dbl> 30.11938, 19.00391, 25.03414, 23.…
## $ Health.Antibiotic.3months <lgl> FALSE, TRUE, TRUE, TRUE, FALSE, F…
## $ Sexual.Intercourse.24hours <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ plate <chr> "PlateBA", NA, "PlateBB", "PlateB…
## $ volume <dbl> 2.2025248, NA, 1.5580467, 1.83601…
## $ dna_conc <dbl> 11.169, NA, 15.789, 20.648, 24.65…
## $ read_conc <dbl> 17233.4038, NA, 20634.8110, 4197.…
## $ read_conc_norm <dbl> 0.9693078, NA, 1.1606229, 0.57860…
## $ lib_size <dbl> 37957, NA, 32150, 7707, 24123, 29…
## $ `Sample month` <dbl> 7, NA, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
## $ `Transit time` <dbl> 1, NA, 1, 1, 1, 2, 1, 1, 6, 2, 1,…
## $ valencia_subcst <chr> "I-A", "IV-C0", "III-A", "IV-B", …
## $ valencia_cst <chr> "I", "IV-C", "III", "IV-B", "IV-B…
## $ max_genus1 <chr> "L. crispatus", "Prevotella", "L.…
## $ max_genus2 <chr> "Gardnerella", "Other", "Prevotel…Continuamos con la tabla “taxa”:
glimpse(vdata$taxa)## Rows: 942
## Columns: 12
## $ taxon <chr> "t1005", "t1019", "t1023", "t1025", "t1030", "t1031"…
## $ taxon_id <chr> "t1005", "t1019", "t1023", "t1025", "t1030", "t1031"…
## $ kingdom <chr> "Bacteria", "Bacteria", "Bacteria", "Bacteria", "Bac…
## $ phylum <chr> "Bacteroidetes", "Bacteroidetes", "Bacteroidetes", "…
## $ class <chr> "Bacteroidia", "Bacteroidia", "Bacteroidia", "Bacter…
## $ order <chr> "Bacteroidales", "Bacteroidales", "Bacteroidales", "…
## $ family <chr> "EU845084_f", "Porphyromonadaceae", "EU845084_f", "E…
## $ genus <chr> NA, "Porphyromonas", "DQ003626_g", NA, "Parabacteroi…
## $ species <chr> NA, NA, NA, NA, "johnsonii", "merdae", NA, NA, "putr…
## $ sequence <chr> "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAG…
## $ genus_oldtaxonomy <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ genus_nosubgroups <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …Ahora veamos la tabla “counts”:
glimpse(vdata$counts)## Rows: 6,714
## Columns: 3
## $ count <dbl> 15, 80, 8, 3, 7, 6, 20, 72, 49, 65, 4, 45, 5, 2, 18, 140, 2,…
## $ sample_id <chr> "s470", "s879", "s894", "s1145", "s1218", "s1388", "s1631", …
## $ taxon_id <chr> "t1005", "t1005", "t1005", "t1005", "t1005", "t1005", "t1005…Luego echamos un vistazo rápido al número total de muestras (n_samples), ASV (n_taxa), y lecturas (n_reads) en el objeto tidytacos:
tacosum(vdata)## n_samples n_taxa n_reads
## 200 942 5022096Hagamos un diagrama de barras apilado de un subconjunto de muestras
Podemos crear muy fácilmente un gráfico para explorar un subconjunto de nuestras muestras (por ejemplo, solo muestras de participantes que han tomado antibióticos en los últimos 3 meses y son menores de 40 años) de la siguiente manera:
vdata %>%
filter_samples(Health.Antibiotic.3months == TRUE, General.Age <= 40) %>%
tacoplot_stack()
La función filter_samples hace lo que dice: filtrar muestras. También
eliminará los taxones de la tabla de taxones que tengan cero lecturas
totales en las muestras restantes. La función tacoplot_stack devuelve
una buena visualización de diagramas de barras apiladas de los taxones
más abundantes en nuestras muestras.
¿Se te ocurren otras ideas de como filtrar y rápidamente visualizar vdata?
Subconjunto de datos
Nuestra siguiente pregunta para este conjunto de datos es hasta qué punto la actividad sexual dentro de las 24 horas antes de tomada las muestras afecta la composición microbiana de las participantes en edad reproductiva (asumamos <=40 años). Para hacernos una idea primero lo podemos visualizar:
vdata_less40 <- vdata %>% filter_samples(General.Age <= 40)
tacoplot_stack(vdata_less40)+
geom_point(aes(y=-0.02,color=Sexual.Intercourse.24hours)) +
geom_point(aes(y=-0.05,color=valencia_cst))
## Diversidad Alpha Para explorar la diversidad alfa, creemos una
versión enrarecida (rarefied en inglés) del conjunto de datos:
vdata_rar <- vdata %>%
add_total_count() %>%
filter_samples(total_count >= 2000) %>%
rarefy(2000) %>%
add_alpha()La función add_total_count agregará números totales de lectura de muestra a la tabla de muestra.
La función rarefy submuestreará aleatoriamente todas las muestras n veces. Solo funciona si el recuento de lecturas de cada muestra es igual o superior a n.
Para determinar la riqueza de ASV, optamos por enrarecer primero, pero esto puede depender de sus datos.
La función add_alpha se puede utilizar para agregar varias métricas de diversidad alfa a la tabla de muestra.
Ahora recalculemos el conteo total de lecturas por cada muestra. Debería ser de 2000 lecturas para todas las muestras.
vdata_rar %>%
mutate_samples(old_total_count = total_count) %>%
select_samples(-total_count)%>%
add_total_count()## $counts
## # A tibble: 5,055 × 3
## count sample_id taxon_id
## <int> <chr> <chr>
## 1 2 s470 t1005
## 2 12 s879 t1005
## 3 4 s1631 t1005
## 4 7 s1817 t1005
## 5 2 s1902 t1005
## 6 4 s1942 t1005
## 7 4 s2790 t1005
## 8 1 s100 t1019
## 9 1 s1902 t1023
## 10 10 s1331 t1025
## # ℹ 5,045 more rows
##
## $samples
## # A tibble: 200 × 32
## sample sample_id Technical.run_20201009_isala…¹ Technical.run_202010…²
## <chr> <chr> <lgl> <lgl>
## 1 SAMEA13411869 s16 FALSE FALSE
## 2 SAMEA13410951 s100 TRUE FALSE
## 3 SAMEA13411960 s114 FALSE FALSE
## 4 SAMEA13412020 s170 FALSE FALSE
## 5 SAMEA13410958 s177 TRUE FALSE
## 6 SAMEA13412032 s183 FALSE FALSE
## 7 SAMEA13412067 s222 FALSE FALSE
## 8 SAMEA13412068 s223 FALSE FALSE
## 9 SAMEA13412130 s288 FALSE FALSE
## 10 SAMEA13412143 s302 FALSE FALSE
## # ℹ 190 more rows
## # ℹ abbreviated names: ¹Technical.run_20201009_isala_cross_1,
## # ²Technical.run_20201016_isala_cross_2
## # ℹ 28 more variables: Technical.run_20201023_isala_cross_3 <lgl>,
## # Technical.run_20201030_isala_cross_4 <lgl>,
## # Technical.run_20201106_isala_cross_5 <lgl>,
## # Technical.run_20201120_isala_cross_6 <lgl>, …
##
## $taxa
## # A tibble: 688 × 12
## taxon taxon_id kingdom phylum class order family genus species sequence
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 t1005 t1005 Bacteria Bacteroide… Bact… Bact… EU845… <NA> <NA> TACGGAG…
## 2 t1019 t1019 Bacteria Bacteroide… Bact… Bact… Porph… Porp… <NA> TACGGAG…
## 3 t1023 t1023 Bacteria Bacteroide… Bact… Bact… EU845… DQ00… <NA> TACGGAG…
## 4 t1025 t1025 Bacteria Bacteroide… Bact… Bact… EU845… <NA> <NA> TACGGAG…
## 5 t1030 t1030 Bacteria Bacteroide… Bact… Bact… Porph… Para… johnso… TACGGAG…
## 6 t1040 t1040 Bacteria Bacteroide… Bact… Bact… Riken… Alis… <NA> TACGGAG…
## 7 t1045 t1045 Bacteria Bacteroide… Bact… Bact… Riken… Alis… putred… TACGGAG…
## 8 t1049 t1049 Bacteria Bacteroide… <NA> <NA> <NA> <NA> <NA> TACGGAG…
## 9 t1089 t1089 Bacteria Proteobact… Delt… Desu… Desul… Bilo… wadswo… TACGGAG…
## 10 t1090 t1090 Bacteria Proteobact… Delt… Desu… Desul… Bilo… <NA> TACGGAG…
## # ℹ 678 more rows
## # ℹ 2 more variables: genus_oldtaxonomy <chr>, genus_nosubgroups <chr>
##
## attr(,"class")
## [1] "tidytacos"Podemos visualizar la diversidad alpha de muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no:
vdata_rar %>%
samples() %>%
ggplot(aes(x = Sexual.Intercourse.24hours, y = observed, fill = Sexual.Intercourse.24hours)) +
geom_boxplot(outlier.shape = NA)+
geom_jitter(height = NULL)
Análisis de coordenadas principales (PCoA)
Nos gustaría abordar las diferencias entre muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no. También estamos más interesados en los géneros que en los ASV.
Una PCoA podría ofrecer información:
vdata_genus <- vdata %>%
aggregate_taxa(rank = "genus")
tacoplot_ord_ly(vdata_rar, Sexual.Intercourse.24hours, samplenames = sample,
dim = 2)## Warning in as.dist.default(xdis): non-square matrixLa función aggregate_taxa fusiona todas las filas de la tabla de taxones en un nivel taxonómico específico, en este caso el nivel de género. Como ocurre con todas las funciones de tidytacos, todas las demás tablas del objeto tidytacos se ajustan en consecuencia.
La función tacoplot_ord_ly determinará la abundancia relativa de taxones en las muestras y luego utilizará las diferencias de Bray-Curtis para ordenar muestras en un espacio bidimensional (o tridimensional) según su composición taxonómica. La adición argumental de plotly “_ly” hace que la figura sea interactiva, lo cual es realmente bueno para el trabajo exploratorio. Esto también funciona para otras funciones de visualización.
Relación entre la composición de la comunidad microbiana y variables.
La siguiente pregunta lógica es hasta qué punto la actividad sexual reciente (dentro de las 24 horas) determina la variabilidad de la composición de la comunidad microbiana. No olvidemos que la composición microbiana vaginal varia con la edad, por eso incluyamos edad en el modelo.
perform_adonis(vdata_genus, c("General.Age", "Sexual.Intercourse.24hours"))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = as.formula(paste("counts_matrix", formula_RHS, sep = " ~ ")), data = metadata, permutations = permutations)
## Df SumOfSqs R2 F Pr(>F)
## General.Age 1 1.717 0.02676 5.5016 0.001 ***
## Sexual.Intercourse.24hours 1 0.956 0.01489 3.0607 0.015 *
## Residual 197 61.500 0.95835
## Total 199 64.173 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La función perform_adonis realizará un ANOVA
PERMutacional para determinar el efecto de las variables de ¨samples¨ en
las diferencias de Bray-Curtis de las comunidades. El resultado muestra
que la edad es un contribuyente a la composición de la comunidad
microbiana (R cuadrado = 0.02676).
Análisis de abundancia diferencial
A continuación, nos gustaría saber cuáles de los 20 géneros más abundantes son significativamente más abundantes en las participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no.
vdata_genus <- vdata_genus %>% add_codifab(Sexual.Intercourse.24hours,
max_taxa = 20)
vdata_genus$taxon_pairs <- filter(vdata_genus$taxon_pairs, wilcox_p < 0.05)
tacoplot_codifab(vdata_genus, FALSE_vs_TRUE)
La función add_codifab agregará una tabla llamada
taxon_pairs al objeto tidytacos, con la abundancia diferencial del taxón
entre las dos condiciones (con respecto al taxón de referencia), para
cada par de un taxón y un taxón de referencia.
La función tacoplot_codifab devuelve un gráfico para
visualizar la abundancia diferencial de taxones entre condiciones, en
comparación con todos los demás taxones como referencia. Podemos
observar que es más probable que Staphiloccus y L. iners group sea
típico encontrar cuando una participante ha tenido actividad sexual
dentro de las 24 horas antes de tomada las muestras.
Es de destacar que existen muchos métodos de análisis de abundancia diferencial y ninguno de ellos es perfecto. Interprete tus resultados con cuidado.
Descargo de responsabilidad: Esta guía es una adaptación de la versión en inglés del tutorial de tidytacos.